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Model equations describing large scale buoyant conveclion in an enclosure are formulated with the vorticity 
and stream function as dependent variahles. The mathematical model, based on earlier work of the authors, is 
unique in two respects. First, it neglects viscous and thermal conductivity effecls. Second the fluid is taken to 
be thermally expandable: large density variations are allowed while acoustic waves are filtered out. A volumetric 
heal source of specified spatial and temporal variation drives the flow in a two-dimensional rectangular enclosure. 
An algorithm for solution of the equations in this vorticity, stream-function formulation is presented. Results of 
compulations using this algorithm are presented. Comparison of these results with those obtained earlier by the 
authors using a finite difference code to integrate the primitive equations show excellent agreement. A method 
for periodically smoothing the computational results during a ealeulalion, using Lanczos smoothing, is also 
presented. Compulations with smoothing al different lime intervals are presented and diseussed. 
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1 . Introduction 

Over the past few years, the National Bureau of Standards has sponsored a joint researeh projeet between 
the Center for Fire Research and the Center for Applied Mathematics to develop, starting from basic con- 
servation laws, a mathematical model of fire development within a room. Large scale convection is an essential 
component of such a model because this fluid motion governs the smoke and hot gas transport within a room 
and also supplies fresh oxygen to the fuel to sustain eombustion. Therefore, development of a mathematical 
model of buoyant convection was begun as a first step towaid a more complete room-fire model, which would 
include effects of combustion chemistry, radiation, and smoke dynamics. The mathematical model for con- 
vection, the partial differentia] equations, and boundary conditions, arc derived in reference [l]. 1 

As noted in earlier papers [1,2] the mathematical model is unique in two respects. First, it is assumed 
that viscous and thermal conductivity effects are negligible. Second, the fluid has been taken to be thermally 
expandable so that large temperature and density variations can be taken into account, while acoustic waves 
have been filtered out to reduce computational time. 

The model equations were integrated for density, pressure, and velocity components by finite difference 
techniques; the algorithm is presented in detail in reference [2]. The algorithm has been verified by comparison 
with solutions to the equations in special cases obtained by analytical and independent numerical means; 
the verification is described in references [3] and [4] and in the present study. 

In section 2 the model equations are recast into a form such that the dependent variables are density, 
pressure, vorticity, velocity potential, and stream function. This formulation, the so-called vorticity, stream- 
function formulation, is an alternate one to that described in reference [2], which we call a primitive-variables 
formulation. An algorithm for integration of the equations in a vorticity, stream-function formulation is also 
presented in this section. 
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Results from the two algorithms are compared: the results are in such good agreement that the difference 
cannot be seen on plots of the dependent variables. This comparison represents a final check on the validity 
of the integration algorithms and their computer code implementations. Therefore, we believe that the solutions 
obtained by these computations are excellent approximate solutions to the model partial differential equations 
presented in reference [l]. 

The model has been developed for two-dimensional, time-dependent fires evolving in a room (rectangular 
enclosure). The fire has been modeled as a volumetric heat source of specified spatial extent and temporal 
behavior. In section 3 density or temperature (which are inversely related at any specified time in this model), 
vorticity, pressure, and velocity potential plots at fixed times during the heating are shown for a sample 
computation. 

It is well known that, because of the quadratic nonlinearity in convection, initially smooth flow fields 
become increasingly more fnrrowed as time progresses; i.e., energy cascades from lower wavenumber modes 
to higher ones. The computational results display this behavior, and the flow field becomes more intricate 
with increasing time, the resolution of the grid providing the only limitation to the resolvable detail. However, 
such an accumulation of energy at a wavenumber inversely proportional to the grid size is both unphysical, 
and, if the computation is carried out long enough, disastrous. Local gradients of the dependent variables 
become much too large and the computation ultimately fails. Therefore, also shown in section 3 of this paper 
are preliminary considerations of smoothing or filtering of the computational results. Such smoothing acts 
analogously to viscosity and can be used to prolong the lifetime of the computation. A brief discussion of the 
effects of a particular type (Lanczos) of smoothing is presented, and resnlts obtained using this smoothing 
are shown. 

2. Formulation 

2.1. Continuous equations 

In an earlier paper [1], the authors had derived a set of nonlinear equations describing very nonadiabatic 
buoyant flows of a nondissipative perfect gas. The magnitude and the spatial variation of the heat source 
(representing the exothermic reaction in a fire) were taken as known. The fluid and the fire source were 
assumed confined in a closed rectangular room with the center of the source along the floor. In contrast to 
reference [l], in this paper we consider only a completely enclosed room (no leaks), and when difference 
equations are introduced, we confine attention to the two dimensional evolution of the flow. 

In this section, the equations derived in reference [l] for a two-dimensional configuration are rewritten so 
that vorticity and the stream function are primary variables, and the finite difference methods used to solve 
the revised equations are presented. 

Equations (11) of reference [l] are 
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Here p is density, u ; the velocity in the i' h coordinate direction (i = 1, 2, 3), p is the pressure excess above 
the mean pressure p a (t) in the room, T the temperature, C p the constant-pressure specific heat, R the gas 
constant, k t g is the gravitational acceleration, and Qipi^i) the specified volumetric heat source. The spatially 
uniform mean pressure p (t) depends only upon time and increases beeause of the heating within the room. 
It is determined in a completely enclosed room by the equation 
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where 7 is the ratio of specific heats, V is the volume of the room and the integration is performed over this 
entire volume. 

We take the substantial derivative of the equation of state and use this with the energy equation to eliminate 
the temperature. The resulting equation describes the evolution of the density under heating 
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Equation (3) and the continuity equation identify D(x i ,t) as the divergence 






(5) 



Finally, as in reference [1], the equation for the spatially variable portion of the pressure is obtained by 
dividing the momentum equations by density and taking the divergence of these equations. The resulting 
equation is 
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The boundary conditions on these equations are that velocity normal to any (impermeable) wall vanish. 

U& = (7) 

where n t are the normal components of a vector describing the boundary wall. From eqs (1) and these 
conditions, the appropriate boundary conditions on the pressure equation are obtained 
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In two dimensions (no dependence on 2), these equations become 
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where x and y are the horizontal and vertical coordinates with velocity components u and v respectively and 
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dx dy 
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is the only nonzero component of vorticity. 

Equations (9) can be recast into a form in which the vorticity, the stream function and the velocity potential, 
together with density and pressure, are primary dependent variables. The velocity components can be written 
as 



u = — + — 
Bx dy 

dy dx 



(11) 



where <]) is the velocity potential and 4 1 the stream function. Equations for <j) and (Jj come from the divergence 
eq (5) and from the definition for vorticity (10) 



V 2 <]> = D(x,y,t) 
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For these two elliptic equations, the stream function and the normal derivative of the potential are zero on 
the boundary: 
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The equation describing the density evolution remains as it was in eqs (9), and an equation for the vortieity 
evolution comes from taking the curl of the two velocity equations in eqs (9). 
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In figure 1 a schematic diagram of a fire evolving in a room in two dimensions and a set of coordinate axes 
are shown. It is assumed that initially the enclosure is filled with quiescent, stratified fluid of density p„(y). 

Smoke and 
hot gases 




*- x 



FIGURE 1. Schematic diagram of a Iwo dimensional fire evolving in a room: 
il is assnmed that ihere is no dependence upon z of any of the properties of 
the fire or of the induced flow field. The fire, ioealized along the floor, has 
a plume of combustion products rising above il. The smoke and hot gases 
rise to the eeiling and fill the room from the top down. 
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We define a density difference from ambient and a pressure difference as follows 

p(*.r,0 = p(x,y,t) - p D (y) 

p{x,y,t) = p(x,y,t) - p£t) + gSlPa(y') dy 

These differences p and p need not be small compared with pjy) and pjt) respectively. Then the first and 
last of eqs (9) and eq (14) become 
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Equations (11), (12) and (15) constitute the complete set of differeutial equations for numerical integration 
in the vorticity, stream function formulation. The boundary conditions are given by eqs (8), with p and p 
replacing p and p, and by eqs (13), 

Finally, we form dimensionless equations using the density p OD = p o (0), the height of the room H as the 
length scale and the free fall time (H/g) l/2 as the time scale. Then, denoting dimensionless quantities with a 
hat 
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Equations (11), (12) and (15) remain exactly the same in dimensionless form with g set equal to one. 
Subsequently, in this paper all quantities will be understood to be dimensionless, and the hat notation will 
be dropped. For the dimensionless coordinates, we note that =£ x =£ \IAR and ^ y =£ 1 where AR = 
H/L. 

2.2. Discrete Equations 

2.2.1. The Basic Algorithm 

In this section the finite difference equations and the boundary relations for the solution algorithm are 
presented. In figure 2a, the two-dimensional rectangular enclosure in dimensionless variables is shown together 
with a schematic representation of the spatial grids used for the finite difference scheme. The grid formed 
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FlGl'fiE 2a. Rectangular enclosure in dimensi'oniess variables ^ * =S , 

AR 

< y ^ 1. The mesli upon which ihc difference seheme is based is shown 

schematically for (/ = J = 4) as a grid of solid lines. The mesh of dashed lines 

joins the center points of the basie mesh cells and is the grid upon which the 

pressure compulation is performed. 



from solid lines represents the basic mesh into which the enclosure is divided: in general there are / mesh 
cells in the k- direction and J mesh cells in the y-direetion. 

Upon this basic mesh, the two components of the vector (u,v) and single component of the vector vorticity 

(0 = — are denned. 

dx dy 

The second grid, formed by joining the center points of the basic grid cells and denoted by dashed lines, 
is that upon which scalar quantities such as density p and pressure p are defined. In Figure 2a the densities 
in the left-hand column of cells and in the bottom row of cells are shown to indicate how they are enumerated 
for the numerical computation. 

In figure 2b a typical mesh cell is shown, illustrating where all of the dependent variables in the finite 
difference scheme are defined relative to the cell. 
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FlCUHE 2b. A typical mesh cell, wilh center tocmed al * = (i — 1/2) S* 
and y = {j — 1/2) Sy, illustrating where all dependent variables for (tie 
finile difference scheme are defined. 
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The following diseretely evaluated functions will denote approximations to the corresponding solutions to 
eqs (11), (12) and (15): 

<j>£ = (H(;-l/2)5*,(/'-l/2)5 r ,n80 

u^= u(ibx,(j-l/2)by,n&t) 

^= ((£- 1/2)8* JSy.nSt) (17) 

p£ = p((i- 1/2)5*, (/-l/2)6y,n8t) 

Pl = p{(i-ll2)hx,{j-U2)hy,nht) 
DJ = D((i- 1/2)6^,0'- l/2)6 r ,n8() , 
cojj = <o(i8:tj8y,nSt) 



where 8* = l/(/-AR) aud Sy = 1// are the mesh cell sizes in the x- and y-directions respectively and where 
8t is the time-step size. Such a staggered grid is commonly used for multidimensional finite difference 
integrations [5], 

With this notation, the following set of finite difference equations was used to approximate eqs (11), (12) 
and (15) and boundary conditions (8) and (13): 

For the first of eqs (15), 1 ^ i ^ /, 1 =S j sS / and ft 5* 2, 



+ (1/2JDSP.G-))} 



sfhere 



P?j ~ P?j ~~ Po(j) = tne density difference from the initial density, 

pjj) = exp[ — (/'— l/2)SyKj = the ambient, initial stratification, (19) 

Y s = the stratification length scale. 

The flux terms F^, and FJ^. for 1 =S i s= /, 1 sS ; =S / are given by 

F „ = PoO'+i) + p.0-1) + P",+i + P",-i fa ~ K,-i 

, p.O'+i) - p.0--i) + p?, + i - ru-i fa + «g,-i N j (20b) 
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For the second of eqs (15), 1 =£ i =S /, 1 =S j ^ J and n Ss 2, 






(21) 



where 
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Equation (18) employs a second-order accurate temporal discretization which eliminates instability that would 
arise if leap frog had been applied. The first of eqs (15) has an undifferentiated term pD{x,y,t) that is well 
known to lead to a computational instability for ordinary differential equations when leap frog is used. Reference 
[10] discusses the simple change in temporal differeneing used here to eliminate this instability. 

Equation (21) uses a straightforward leapfrog temporal differeneing, and both eqs (18) and (21) are started 
by using the same spatial discretization and an explicit, first-order time step. 

At eaeh time step, after the vorticity has been updated, three elliptic equations must be solved, eqs (12a), 
(12b) and the last of eqs (15). Equations (12a) and (12b) are differenced using a standard Five point star 

^i Oft+ij - 2*s + #-!,) + ^5 (4>tj + i - 2fy + 4>?,-i) = D?j (24a) 

— (W +1J - 2^ + ^- lj ) + ^ (K + i - 2 ^ + <-i) = - < (24b) 

and the boundary conditions (13) are introdueed in the usual fashion. These equations have been solved 
using software routines from FISHPAK [6]: eq (24a) was solved using BLKTRI and more reeently using 
POISTG, while eq (24b) was solved using PWSCRT and more recently GENBUN. Routines BLKTRI and 
PWSCRT have limitations on the number of mesh points or unknowns which they can solve, whereas POISTG 
and GENBUN were produced more reeently and are free of sueh limitations. Most computations were performed 
with the former routines, but recently several computations have been performed using the latter. 

The velocities are then obtained from the potential and stream function by difference forms of eqs (11): 



£<+T + i, -«*) + £■ 



u 5 = ^ (*?+u - *s) + r ^ ~ <--i) ( 25a) 



«S = - (4>& +I - 4>3) - - (iP5 - *?_ w ) (25b) 
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For the third of eqs (15), for 1 « i < / and 1 « j « J t 
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where the fluxes F°j and F^j are defined as follows: 



for 1 ^ i ^1 - i,i^j«s; 
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and for 1 *Z j < 7 - l,Ui«/ 
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(27b) 



and where 
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Note that boundary conditions (7) on the normal velocities imply that u j = U/j = f or 1 =S j ^ / and 
v i0 = v t j = for 1 ^ i ^ /. These boundary conditions are applied formally in the expressions for the 
fluxes F^.., Fjjj, , F". and F".. in mesh cells adjacent to boundaries. The boundary conditions (8) in discrete 
form become 



PSj = Pi, 



for 1 < j i =£ / 



(28a) 



Ph = P? +lj 
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Pli - Plo = 8 r(P"i + P"oV2 

for 1 *£ j « / (28b) 

Pu + i ~ Pu = ~ MP"; + i + PI-:,V2 

We note that eqs (26) together with boundary conditions (28) constitute a singular linear algebraic system of 
equations. When eqs (26), with boundary conditions incorporated, are summed, the left hand side sums to 
zero, demonstrating that all of the equations are not linearly independent. Also, the last three terras on the 
right hand side sum to sero, producing the requirement that the double sum over (Z?jJ +1 — D"~ x )l2ht must 
vanish. Examination of eq (30) below for D^ shows that it has been chosen so that its double sum oyer all 
mesh points vanishes, and that the eondition which must be satisfied to allow this choice produces eq (31) 
below for the mean pressure. Therefore, the singular linear algebraie system is seen to be consistent and 
thus to allow a solution. The solution is unique by noting that the double sum over all mesh points of p"j must 
be zero. 

At eaeh time step it is necessary to ealeulate the solution of the linear algebraic system for the pressure, 
eqs (26) with boundary conditions (28) incorporated. The method of solution must take into account non- 
uniqueness of the algebraic system. The solution method must also be able to solve large linear systems 
accurately, sinee there are 1J equations and cumulative errors from many time steps may destroy the com- 
putation. Finally, it is very important that the solution be obtained quiekly since the calculation is made at 
each time step, and hundreds of time steps must be taken. 

The solution method we have adopted is a hybrid method which eombines an iterative algorithm, conjugate 
gradients, with a fast dircet Poisson solver. The eonjugate gradients algorithm provides an iterative technique 
for solving the linear algebraic system of equations. Details of the algorithm are presented in reference [7]. 

The heat source has been chosen to have the form 



Q tj = (^J \ exp [-p(* ( -*«) a - \yj] (29b) 



x t = (i- 1/2)8* , yj = (/-l/2)Sy (29c) 

P = Q a tanh At* (29d) 

B-l 

i° = , t" = 2 8t "' ( 29e ) 



Henee, the diserete divergence of the velocity field becomes 



^■ = 7;[(7-i)<?v-Kir < 30a ) 

tro 

* = ^ £ 2 &■ (30b) 

JJ i = 1 -J = 1 
and the mean background pressure is found from the difference equation 

p n + l = p2 -l + Kfn2htr (31) 

wifhp° = pi = 1 since/' = 0. 
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The linear stability of the algorithm is the only other consideration for diseussion. A linear stability analysis 
of eq (18) for the density shows that the time step 8 J must satisfy the following condition for stability 



8(" 



max 
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m z 
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)u$\ [vgj_ 



-1/2 



(32) 



where 






(33a) 
(33b) 



When the stability condition, eq (32), is not satisfied by a time step, the time step 8i" is halved. Then the 
time-marching algorithm is restarted using the last time-level values as initial conditions. A first-order time 
step is taken and then leap-frog is resumed. 

A linear stability analysis of the difference equation for the vorticity, eq (21), yields exactly the same form 
for the stability criterion as that found above for the density equation. Reference to figure 2b shows that the 
density and vorticity are evaluated at different points in the mesh, however, and therefore, the divergence D 
and the velocity components U and V are to be evaluated at different points than those used in eq (32), To 
account for the difference in the stability criterion implied above, in all calculations performed using the 
algorithm described above, the time step was chosen to be less than or equal to 0.8 the maximum value 
found for the right hand side of eq (32). 



2.2.2. Lancxos Smoothing 



The nonlinear nature of the equations of fluid dynamics implies that initially smooth data will, in general, 
produee flow fields with fine structure. Since the results presented are for finite difference computations, the 
resolution of the flow field is limited by the grid size used to perform the computation: structures of a size 
comparable to a few grid cells can be resolved, whereas smaller structures may represent artifacts of the 
calculations. In addition, in the computations it has been found that the calculation will eventually fail because 
of the intricate detail (and sharp gradients this detail represents) if the number of time steps becomes large 
enough. 

It is for these reasons that various methods of smoothing data generated by the computations have been 
examined. In the method presented here, the computation is stopped periodically, with a period specified as 
input to the computation, and the data are smoothed spatially. The computation is restarted using the smoothed 
data as initial conditions, the results not being allowed to develop the intricate detail it might otherwise 
develop. The method used to smooth the data is a variation of one suggested by Lanczos in reference [8], 
Using this method with a relatively long smoothing period, computations Wave been extended indefinitely. 

The smoothing used here is that proposed by Lanczos, but is modified slightly for our purposes. In reference 
[8], the smoothed data is obtained from the value of the data at each point by adding a specified multiple 
of the fourth difference at that point. The change in value between the smoothed and unsmoothed data then 
is of order A 4 where h = \IJ and J is the number of mesh points in one direction in space. 

Since the computational scheme described here is only second order accurate in the spatial mesh size, 
0(/i 2 ), a less refined smoothing was used. The smoothing is accomplished by adding a specific multiple of 
the second-order difference at the point to the value of the datum at that point and is 0(/i 2 ). When the method 
is generalized to two dimensions, it becomes equivalent to adding one fifth of the finite difference, five-point 
Laplacian to the value of the datum at each point to obtain the smoothed datum. (This is also equivalent to 
replacing the value at a point by the average of its value and the values of its four nearest neighbors.) 

Therefore, after a specified number of time steps m, the density and vorticity data at time level m are 
changed according to the following prescription: 
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P? <- P* + 5 (pr+u + P?-i., + PS- + i + PO-i - 4p*) (34a) 

»* «-»£ + g «+i,y + wr-ij + *>£ +1 + «£-i - 4cojO (34b) 

At the boundaries, the following rules are used for the smoothing. For cells adjacent to boundaries, the density 
is assumed to have the same value in a ficticious cell outside the boundary as its value in the cell under 
consideration. (This is the difference equivalent of saying that the normal derivative of the density at a 
boundary is zero.) The vorticity is taken to be zero on the boundary, and represents a free slip boundary 
condition, 

A rough estimate of an equivalent Reynolds number or Grashof number introduced by smoothing can be 
made by the following argument. The effect of viscosity in the vorticity equation in the Boussinesq approx- 
imation arises, when the equations are made dimensionless in an appropriate way, as — - V 2 co, where Re is 

Re 

the Reynolds number, V 2 is the Laplacian and to is the vorticity. The average effect per time step of Lanczos 

ch? 
smoothing in the vorticity equation can also be represented as — r V 2 <t> where c is a constant of order unity, 

mo 

5 is the time step, h is the mesh spacing, m is the number of steps between smoothings and V 2 is the 

discretized, five-point representation of the Laplacian. Equating the coefficients nf the Laplacian operators 

provides an estimate of the Reynolds number introduced by smoothing: 

8 m 8 

Re = -— = -ml -J 

c hr c 

where / and / are the number of mesh points in each direction. Taking c = 1, 5 = 0.1, m = 40, / ■ / s 
1000, then Re = 4 X 10 3 , and noting that the Grashof number Gr is approximately the Reynolds number 
squared, Gr = 1.6 X 10 7 . 

3. Computational Results 

As discussed in the Introduction, the vorticity, stream-function algorithm and a code implementing this 
algorithm were developed as a method for solving the partial differential equations derived in reference [l]. 
The other method for solving these equations, a finite difference method for directly integrating the equations 
of motion in primitive variables (density, pressure and the two components of velocity) was described in 
reference |2]. Reference [3] describes comparisons of the results computed using the primitive variable 
algorithm with analytical results obtained in special cases. These comparisons were performed to test the 
algorithm and the computer-code implementation. Final comparisons were made between results computed 
using the primitive variables code and those eomputed using the vorticity, stream-function code. Agreement 
between results was found to about five significant figures after a few time steps and to between three and 
four significant Figures after hundreds of time steps. The discrepancy between results is smaller than the 
errors introduced by discretization for the mesh sizes used and well below differences which eould be observed 
by plotting. 

In this section some computational results are presented and discussed. The density, pressure and vorticity 
are sealar functions of the horizontal and vertical coordinates at any specified time. We have found that 
contours of eonstant value of any of these scalar quantities are a convenient way to display them. Since the 
temperature and density are inversely related at any particular time, eontours of eonstant density are also 
contours of eonstant temperature. 

All contour plots were prepared from a graphics package developed by the National Center for Atmospheric 
Researeh. The numbers indicating contour values are relative only. Solid lines represent values of the variable 
greater than zero and dashed lines represent values less than zero. For the results presented here, density 
and temperature are inversely related: eontours of constant density have been labeled as contours of constant 
temperatures for illustration. Therefore, for the temperature plots, temperature contours above a reference 
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value appear dashed and below that reference are solid. Graphic details on the scale of a mesh cell, which 
is determined by the distance between fiducial marks on the sides of the plots, should be ignored. 

All computations have been performed on one of three computers, the NBS UN1VAC 1108, the U.S. 
Treasury UN1VAG 1100/81 or the Cybernet Cyber 175. The computations require about 45K words of storage 
for the 31 X 31 mesh and were performed on any one of the three machines. Typical running times on the 
1108, which is the machine most frequently used, is about 30 to 45 min of CPU time for 200 to 300 time 
steps. 

In figure 3 contours of constant temperature (isotherms) are shown at four dimensionless times for a 
volumetric heat source centered along the floor in a square room. The rate of heat added per unit volume is 
largest along the floor at the center of the room and decreases in a Gaussian fashion with horizontal distance 
from the center and exponentially with height above the floor. The heat source is "turned-on" as a hyperbolic 
tangent with respect to time asymptoting to full strength around t = 1.0. (Note that T denotes time on the 
figure titles.) At the first time, t = 2, the problem is still linear; the flow velocities are sufficiently small 
than convection is unimportant, and the temperature increase in the fluid is directly proportional to the 
volumetric rate of heat added. Therefore, the isotherms are also contours along which the volumetric heat- 
addition rate is constant. (These contours are found to be parabolas.) 

These computations were performed on a spatial mesh having 31 cells in the horizontal and 31 cells in 
the vertical directions; the tick marks along the boundary of the enclosure show the mesh cell spacing. 

At time 8.5, the second frame of figure 3, a buoyant thermal has developed, giving the appearance of a 
mushroom cloud. The buoyant thermal intensifies in strength until the thermal hits the ceiling, as shown in 
the third frame, t = 11.5 and begins to spread. Inside the plume a distinctly periodic structure has begun 
to develop, as can be seen vividly in this frame; here, progressing up the plume along its centerline, one 
finds a local low first, then a periodic sequence of local highs and lows up to about the center of the head 
of the thermal. 

The heated gases spread along the ceiling and fill the room from the top down, as shown in the last frame 
of figure 3 at t = 14.5. This physical behavior is exactly what is observed in room-fire tests and in other 
experimental observations of heating in enclosures. The symmetry about the centerline of the room displayed 
in these computations is some measure of the accuracy with which they were performed: the heat source is 
placed symmetrically, but the computations were performed as if no symmetry existed. 

In figure 4 contours of constant vorticity at various times are displayed. The contours show an anti-symmetry 
about the centerline, as would be expected for the vorticity, and the physical features described for the density 
(or temperature) contours are reflected in the vorticity plots. Because the vorticity represents a higher order 
difference of the dependent variables, these contours display more fine scale (on the mesh scale) features 
than the density contours. Also, later in the computations, "noise" of a nonsymmetric and fine scale begins 
to show up. The vorticity plots are not as "smooth" as those of the density (or temperature). 

In the first frame of figure 4 at time 2,0, early in the heating process, two vortices of equal magnitude and 
opposite sign develop with centers in the regions of steepest gradient of the heat source. Convection has 
begun by dimensionless time t = 8.5, frame 2 of figure 4, and the vortices are pinched together and buoyed 
upward off the floor. The vortices are convected toward the ceiling. Also vorticity of periodically varying 
strength is generated within the source region and the strength of the vorticity increases with distance above 
the source. Finally, frames 3 and 4 of figure 4 show that the periodically varying vorticity trains split when 
the ceiling is encountered and form two large regions of vorticity of opposite sign. 

The pressure as a function of horizontal and vertical coordinates at any specified time can also be displayed. 
In figure 5 contours of constant pressure at four times during the room filling are shown. This pressure actually 
represents only the spatial variation of the pressure and has been normalized at each time so that its mean 
value is zero. This is the quantity which, together with buoyancy, induces the flow. As with the other contour 
plots, solid lines represent contours with values greater than zero and dashed lines indicate values less than 
zero. The pressure plots are seen to be smoother than those of density (temperature). 

Early, the first frame at t = 0.5 of figure 5, the pressure is highest at the source, where heating takes 
place. As convection starts, the high pressure region is lifted off the floor, time t = 2.0: a significant enough 
convective velocity has developed by dimensionless time t = 2.0, that a low pressure region due to a Bermoulli 
effect can be seen at the floor. This low pressure region is associated with the high convective velocities, or 
the vorticity pair shown in figure 4. The low pressure region develops a double minimum, symmetric about 
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FIGURE 3. Isothenns at four dimensionless times, T, during healing by a volumetric 1ieat source centered along tlie Floor in a square room. These four 
frames illustrate (1) early heating widi no convection, (2) a starting pluiiie, (3) the plume impacting ike ceiling, and <4) the room filling from the lop down. 
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FIGURE 4, Contours of conslanl vorticity al four dimensionless limes, T, during healing by a volumetric heat source centered along ihe floor in a square 
room. These four frames correspond lo Lhe ones shown in figure 3. 
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FIGURE 5. Contours of constant pressure al four dimensionless limes, 7", during healing. The calculation shown here is ihe same as lhat displayed in 
figures 3 and 4, bul ihe times are chosen lo illustrate the stages in the development of the pressure field. 
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the room ccntcrlinc as seen at time 8.5 and rises to the ceiling, as shown in figure 5, time 11.5, where a 
strong compression develops at the center of the ceiling. The final pressure diagram displays the double 
minimum, associated with the strong vortex pair and high temperature shown in figures 3 and 4, separating 
and moving toward the walls. 

In figures 6 and 7 contours of constant potential <t> and constant stream funetion i|i are shown. Only one 
plot of the potential is shown because the spatial dependence does not change with time: the potential function 
is separable in space and time. Four frames of stream function are shown in figure 7. The stream function 
is antisymmetrical about the centerline and displays a peak and a valley which slowly rise toward the ceiling. 
The stream function is rather smooth, showing only a slightly wavy behavior at later times. 
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FIGURE 6. Velocity-potential contours al dimensionless lime, T = 1.0 /or 
the calculation displayed in figures 3—5. The potential function ie a function 
of space times a function of lime; therefore, the spatial dependence does 
not change with time. 



The base computation, shown in figures 3—7, was repeated several times with smoothing introduced at 
different numbers of time steps. Figures 8 and 9 compare constant temperature contours at dimensionless 
time 11.5 and 14.5 respectively determined by the base computation, on the upper left, and by three levels 
of smoothing in the other three plots of each figure. It is seen that the fine structure is eliminated by smoothing, 
but large scale features are still retained. In figure 8 four plots of constant temperature contours at approximately 
the same time are compared. The plot in the upper left hand corner is the unsmoothed computation. The 
plot at the upper right is the result of the computation smoothed only once up until that time, after N = 160 
time steps. The next plot, lower left, is from a computation smoothed every 80 time steps, that in the lower 
right is from a computation smoothed every 40 time steps. Figure 9 shows a similar comparison of the effects 
of smoothing, but at approximately t = 14.5. These plots show clearly the loss in detail with increasing 
frequency of smoothing or increased simulated viscosity. They also show that the buoyant plume rise slows 
due to the decreasing gradients. The smoothing and loss of fine-scale detail with increased frequency of 
smoothing are apparent, and, in fact, these results appear much closer qualitatively to results obtained in 
previous studies which integrated the Navier-Stokes equations by finite difference techniques (for example, 
reference [9], figures 4 and 5). 
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FiCUHE 7. Stream-functicu conlours at lour dimensionless limes corresponding to ihe limes displayed in figures 3 and 4. 
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FIGURE 8. Isotherms al approximately the same dimensionless time, T, with Lanczos smoothing applied at different frequencies. The frame in the upper 
left comer is ihe unsmoolhed computation. The frame at the upper right is smoothed every 160 time steps, thai at the lower left every 80 time steps and 
thai at the lower right every 40 lime steps, 
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FIGURE 9. Isotherms al approximately the same dimensionless lime, T, with Lanezos smoothing applied at different freqnencies. The frame in the upper 
left eorner is the unsmoolhed computation. The frame in the upper right is smoothed every 160 time steps, that al the lower left every 80 time steps nnd 
that nt the lower right every 40 time steps. 
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